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Dynamic critical behavior of multi-grid Monte Carlo for 
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We introduce a new and very convenient approach to multi-grid Monte Carlo (MGMC) algorithms for general 
nonlinear <7-models: it is based on embedding an XY model into the given <7-model, and then updating the 
induced XY model using a standard XK-model MGMC code. We study the dynamic critical behavior of this 
algorithm for the two-dimensional O(iV) <7-models with N = 3,4, 8 and for the SU{Z) principal chiral model. We 
find that the dynamic critical exponent z varies systematically between these different asymptotically free models: 
it is approximately 0.70 for 0(3), 0.60 for 0(4), 0.50 for 0(8), and 0.45 for SU(Z). It goes without saying that 
we have no theoretical explanation of this behavior. 



Multi-Grid Monte Carlo (MGMC) [1-8] is a 
collective-mode approach that introduces block 
updates (of fixed shape but variable amplitude) 
on all length scales. The basic ingredients of the 
method are: 

1) Interpolation operator: A rule specifying 
the shape of the block update. The interpolations 
most commonly used are piecewise- constant and 
piecewise-linear . 

2) Cycle control parameter 7: An integer 
number that determines the way in which the dif- 
ferent block sizes are visited. In general, blocks of 
size 2 l are updated f l times per iteration. Thus, 
in the W-cycle (7 = 2) more emphasis is placed 
on large length scales than in the V-cycle (7 = 1). 

3) Basic (smoothing) iterations: The local 
Monte Carlo update that is performed on each 
level. Typically one chooses to use heat-bath up- 
dating if the distribution can be sampled in some 
simple way, and Metropolis otherwise. 

4) Implementation: The computations can 
be implemented either in the recursive multi-grid 
style using explicit coarse-grid fields [9,1-4], or 
in the unigrid style using block updates acting 
directly on the fine-grid fields [10,5-8]. For a d- 
dimensional system of linear size L, the compu- 



tational labor per iteration is 

Work(MG) ~ L d ify<2 d (1) 

for the recursive multi-grid approach 2 , and 

' L d log L for 7 = 1 
L d + io g2l f orT>1 



Work(UG) 



(2) 



for the unigrid approach. Thus, the unigrid im- 
plementation is marginally more expensive for a 
V-cycle, but prohibitively more expensive for a 
W-cycle. 

The efficiency of the method can be analyzed 
rigorously in the case of the Gaussian (free-field) 
model, for which it can be proven [2,11] that crit- 
ical slowing-down is completely eliminated. That 
is, the dynamic critical exponent z (the expo- 
nent with which the autocorrelation time r di- 
verges as the correlation length £ tends to infin- 
ity) is zero. 3 More precisely, the algorithm with 
piecewise-linear interpolation exhibits z = for 
both V-cycle and W-cycle, while the one with 
piecewise-constant interpolation has z = only 
for the W-cycle (the piecewise-constant V-cycle 
has z = 1). 



2 For the W-cycle in d 
logL. 



1 there appears an extra factor 
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3 See [12] for a pedagogical discussion of the various au- 
tocorrelation times and their associated dynamic critical 
exponents. 
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One is therefore motivated to apply MGMC 
to "nearly Gaussian" systems: one might hope 
that critical slowing-down would be completely 
eliminated (possibly modulo a logarithm) or at 
least greatly reduced compared to the z f» 2 
of local algorithms. In particular, we are inter- 
ested in applying MGMC to asymptotically free 
two-dimensional cr-models, such as the TV-vector 
models [also called 0(TV)-invariant cr-models] for 
TV > 2, the SU(N) principal chiral models, and 
the RP N ~ 1 and CP^ -1 cr-models. 

In view of the rigorous results for the Gaussian 
case, we want to investigate two questions: 

1) Is z = for all asymptotically free two- 
dimensional cr-models? If not, does z vary from 
one asymptotically free model to another? 

2) Is the algorithm with piecewise-constant in- 
terpolation and a W-cycle as efficient as the one 
with piecewise-linear interpolation and a V-cycle, 
i.e. is zpc,w equal to zpi, t v for these models? 

The key design choice in a MGMC algorithm 
is that of the interpolation operator; indeed, this 
choice determines most of the remaining ingre- 
dients. If one chooses a "smooth" interpolation 
such as piecewise-linear, then it is usually im- 
possible to implement true recursive MGMC 4 , as 
there is no simple form for the induced coarse-grid 
Hamiltonians. Therefore one is obliged to use the 
unigrid style, a V-cycle, and Metropolis updating. 
We call this the "German" approach [5-8]. On 
the other hand, if one chooses a "crude" interpo- 
lation such as piecewise-constant, then is obliged 
to use a W-cycle in order to have a chance at 
z < 1; but piecewise-constant interpolation usu- 
ally gives rise to a simple coarse-grid Hamiltonian 
(typically a slight generalization of the fine-grid 
Hamiltonian), so that one can use the recursive 
multi-grid style and, at least in principle, heat- 
bath updating. This is the approach taken by 
our group [1-4]. 

The "German" version has the advantage of be- 
ing easy to implement for diverse models, but its 
use of Metropolis updates introduces several free 
parameters that have to be adjusted, making it 
more difficult to test systematically. "Our" ver- 



4 In particular, this is the case for the nonlinear cr-models. 



sion has no free parameters, but its implemen- 
tation is cumbersome and model-dependent, in 
the sense that the program (and in particular 
the heat-bath subroutine) has to be drastically 
rewritten for each distinct model. For example, 
among the TV-vector models, the only ones that 
can be handled conveniently are TV" = 2 [3] and 
TV = 4 [4], by exploiting the isomorphism with 
the U{\) and SU{2) groups, respectively. 

With this problem in mind, we have developed 
a new implementation of MGMC that combines 
some of the advantages of both methods; in par- 
ticular, it can be used conveniently for a large 
class of cr-models with very little modification of 
the program. The idea is to embed angular vari- 
ables {9 X } into the given cr-model, and then up- 
date the resulting induced XY model by our stan- 
dard (piecewise-constant, W-cycle, heat-bath, re- 
cursive) MGMC method. We do not claim that 
this approach is superior in practice to the "Ger- 
man" method — that remains to be determined 
— but we do think that it is well suited for the 
systematic study of the dynamic critical behavior 
of MGMC algorithms. 

For the models discussed here, the induced XY 
Hamiltonian is of the form 

Hembed = — [<X XX > £Os(6 x — 6 x i) + 

("') (3 XX , sm(6 x - 9 X ,)] , (3) 

where the induced couplings {a xx i, f3 xx i} are 
given in terms of the values of the original spins. 
We illustrate this for the two cases thus far stud- 
ied: 

1) TV-vector models: The original variables are 
unit vectors in R^, and the original Hamiltonian 
is 

7~(-N-vector = ~P ^ & x ' & x 1 • (4) 

(sir') 

The embedding is given by choosing randomly a 
plane P in R^, and defining B x to be the angular 
coordinate of the projection of cr x onto P. (At 
each iteration a new random plane is chosen.) If 
each vector cr x is decomposed into its parts paral- 
lel and perpendicular to the plane, a x = a\-\-a x , 
then the induced XY Hamiltonian is of the form 
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(3) with 

<W = (5) 
Pxx> = (6) 

Though the couplings {a xx i} are random, the in- 
duced XY model is ferromagnetic if the original 
TV-vector model is. 

2) SU(N) principal chiral models: The original 
variables are group elements U x £ Sf7(TV), and 
the Hamiltonian is given by 

-h SU {n) = -p Retr (^i u *>) ■ ( 7 ) 

(a; a;') 

In this case we choose to write the embedding 
from a different point of view, namely we define 
the updated variable JJ™ ew by 

fjnew = Re i6 w T R -1 V M ^ (g) 

where R is an SU(N) matrix picked at random 
at each step (the analogue of the random plane 
in the previous case), and T £ Su(TV) is a fixed 
traceless diagonal matrix with entries ±1 or [we 
use T = diag(l, -1, 0, . . ., 0)]. T satisfies 

e i6T = T 2 CQs9 + iT sin6l + (J _ T 2 ) _ ^ 

Therefore, the induced XY Hamiltonian is of the 
form (3), with couplings 

a xx , = /3Retr(*7j RT 2 iT 1 U x .) (10) 
(3 XX , = /3Imtr(*7j RT R- 1 U x .) (11) 

The embedded model is simulated with initial 
condition 6 X = (i.e. U£ ew = U° ld for all x). 
Note that the couplings in this case are not only 
disordered but are in general frustrated. 

We have investigated the dynamic critical be- 
havior of the MGMC method (using piecewise- 
constant interpolation, W-cycle and heat-bath 
updates) for a variety of cr-models in two dimen- 
sions. A few years ago, a careful study [4] of 
the direct MGMC method for the 4-vector model 
showed Zi nt> M 2 = 0-60 ± 0.07 (compared to z ps 2 
for conventional local algorithms such as heat 
bath and Metropolis). Therefore, critical slowing- 
down is greatly reduced but not completely elim- 
inated. Indeed, it is not completely eliminated 



even for the one-dimensional version of the same 
model, for which our data [13] are consistent with 
a logarithmic growth of the autocorrelation time. 

Recently, we have applied the X Y-embedding 
algorithm to the TV-vector model for N = 3, 4, 8, 
finding that: 

1) For TV" = 4, the embedding algorithm be- 
longs to the same dynamic universality class as 
the direct MGMC algorithm, as expected. 

2) The dynamic critical exponent z is TV- 
dependent: we have Zi nt> M 2 ~ 0.70, 0.60, 0.50 for 
TV = 3, 4, 8, respectively (error bars in each case 
are roughly ±0.05). It thus appears that z de- 
creases as TV gets larger; and it might conceiv- 
ably be the case that z tends to zero as TV tends 
to infinity, which would be consistent with the 
vague idea that the TV = oo model is "essentially 
Gaussian" . 

The method was also tested for the 5f7(3) prin- 
cipal chiral model, yielding the preliminary esti- 
mate z = 0.45 ± 0.02. 

These exponents are found by fitting to the dy- 
namic finite- size- scaling Ansatz 

n n t,A{p, L) « i(p, L) Zint - A g A {£(p, L)/L) , (12) 

where A is an observable (here the square of the 
magnetization, which is found to be the slow- 
est mode among those we study), T;„ ()J 4(/3, L) is 
its integrated autocorrelation time [12], £((3,L) is 
the second-moment correlation length described 
in [4], and g A is a smooth function. We plot 
((3, L)/i(P, L) z as a function of £(/?, L)/L, 
and vary z until the points fall nicely on a single 
curve (except for possible corrections to scaling 
for the smaller lattices). Figures 1, 2 and 3 show 
the results for the 0(3), 0(8) and 5f7(3) models, 
using our preferred choices for the exponents z. 
The corrections to scaling appear to be weaker 
for SU(3) than for 0(3) or 0(8); we don't know 
why. 

This work was supported in part by NSF grant 
DMS-9200719. 
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d=2, 0(3) model, z = 0.7 
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Figure 1. T in t,M»{P,I>)/t(J3,L) z versus 

for the 0(3) model, with z = 0.70. Lattice sizes 

L are 32 (+), 64 (x), 128 (□), 256 (O). 



d = 2, 0(8) model, z = 0.o 
10 i i i i I i i i i I i i i i I i i i i I i i i 
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Figure 2. Same, for the 0(8) model, with z = 
0.50. 
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d = 2, SU(3) model, z = 0.45 
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Figure 3. Same, for the 5f7(3) model, with z = 
0.45. Includes also L = 16 (*). 



